///@file Distribution.cpp
///@brief Implementation of the class defining a data distribution
///@author Arnaud Duval
///@version 0.0
///@date 2010/05/12

#include "Distribution.h"
#include "ProbUtils.h"
#include <cstdlib>
#include <cassert>
#include <iostream>
#include <fstream>
#include <ctime>
#include <cmath>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

///@brief Default constructor
Distribution::Distribution()
{
	size = 0;
}

///@brief Constructor with arguments
///@param n size of generated distribution
///@param t type of distribution : "gauss" for a Gaussian distribution, "uniform" for a uniform distrubution"
///@param p1 First parameter for the distribution. If gaussian : avarage value, if uniform : lower bound
///@param p2 Second parameter for the distribution. If gaussian : standard deviation, if uniform : upper bound
Distribution::Distribution(unsigned int n, std::string t, double p1, double p2)
{
	srand((unsigned int)time(NULL));
	if(t == "gauss")
	{
		data = new double[n];
		size = n;
		params = new double[2];
		params[0] = p1;
		params[1] = p2;
		type = t;
		///@brief Here are various methods that must be evaluated
		//Method with conversion of distribution function (not efficient because of the approximation of reciprocal Gauss distribution function)
		//**********************************************************
		/*for(unsigned int i = 0 ; i < size ; i++)
		{
			data[i] = RecFgauss((double)rand()/(double)RAND_MAX, p1, p2);
		}*/
		//**********************************************************

		//Box Muller method, cartesian form

		//Box Muller method polar form
		unsigned int count = 0;
		while (count < n)
		{
			double u1 = (double)rand()/(double)RAND_MAX;
			double u2 = (double)rand()/(double)RAND_MAX;
			double z0 = sqrt(-2.*log(u1))*cos(2.*M_PI*u2);
			double z1 = sqrt(-2.*log(u1))*sin(2.*M_PI*u2);
			data[count] = p2*z0 + p1;
			if(count+1 != n)
				data[count+1] = p2*z1 + p1;
			count += 2;
		}
	}
	else if(t == "uniform")
	{
		assert(p1 <= p2);
		data = new double[n];
		size = n;
		params = new double[2];
		params[0] = p1;
		params[1] = p2;
		type = t;
		for(unsigned int i = 0 ; i < size ; i++)
		{
			data[i] = p1 + ((double)rand()/(double)RAND_MAX)*(p2-p1);
		}
	}
}

///@brief Copy constructor
///@param d Distribution object to copy
Distribution::Distribution(const Distribution& d)
{
	size = d.size;
	data = new double[size];
	for(unsigned int i = 0 ; i < size ; i++)
		data[i] = d.data[i];
	type = d.type;
	params = new double[2];
	params[0] = d.params[0];
	params[1] = d.params[1];	
}
///@brief Copy constructor for temporary objects
///@param d Pointer to Distribution object to copy
Distribution::Distribution(Distribution* d)
{
	size = (*d).size;
	data = new double[size];
	for(unsigned int i = 0 ; i < size ; i++)
		data[i] = (*d).data[i];
	type = (*d).type;
	params = new double[2];
	params[0] = (*d).params[0];
	params[1] = (*d).params[1];	
}

///@brief Destructor
Distribution::~Distribution()
{
	if(size)
	{
		delete data;
		delete params;
	}
}

///@brief Overload of [] operator : direct access to a realization
///@param i Index of returned realization
double& Distribution::operator[](unsigned int i) const
{
	assert((i >=1)&&(i <= size));
	return data[i-1];
}

///@brief Save distribution and its informations to a file
///@param filename Name of the file to write
void Distribution::SaveToFile(std::string filename)
{
	std::ofstream f(filename.c_str());
	f.precision(20);
	f << type << "\n";
	f << params[0] << "\t" << params[1] << "\n";
	f << size << "\n";
	for(unsigned int i = 0 ; i < size ; i++)
		f << data[i] << "\n";
	f.close();
}

///@brief Read distribution and its informations from a file
///@param filename Name of the file to read
void Distribution::ReadFromFile(std::string filename)
{
	std::ifstream f(filename.c_str());
	f >> type;
	params = new double[2];
	f >> params[0] >> params[1];
	f >> size;
	data = new double[size];
	for(unsigned int i = 0 ; i < size ; i++)
		f >> data[i];
	f.close();
}

///@brief Return size of distribution
unsigned int Distribution::Size()
{
	return size;
}

///@brief Return array of data
double *Distribution::Data()
{
	return data;
}

///@brief Return average value of the distribution
double Distribution::Average()
{
	double res = 0.;
	for(unsigned int i = 0 ; i < size ; i++)
		res += data[i];
	return res / size;
}

///@brief Return standard deviation of the distribution
double Distribution::StandardDeviation()
{
	double avg = Average();
	double res = 0.;
	for(unsigned int i = 0 ; i < size ; i++)
		res += pow(data[i] - avg, 2.);
	return sqrt(res/size);
}

///@brief Sort the distribution using bubble sort algorithm
void Distribution::BubbleSort()
{
	assert(size > 0);
	bool permut = true;
	while(permut)
	{
		permut = false;
		for(unsigned int i = 1 ; i < size ; i++)
		{
			if(data[i-1] > data[i])
			{
				double temp = data[i];
				data[i] = data[i-1];
				data[i-1] = temp;
				permut = true;
			}
		}
	}
}

///@brief Return a confidence interval
///@param conf Confidence (in percent)
double* Distribution::Confidence(double conf)
{
	assert(conf >= 0);
	assert(conf <= 100);
	double* res = new double[2];
	conf /= 2.;

	int i1 = (int)((1. - conf/100.)*size);
	int i2 = (int)((conf/100.)*size);
	
	Distribution temp = *this;
	temp.BubbleSort();
	
	res[0] = temp.data[i2];
	res[1] = temp.data[i1];
	return res;
}

///@brief Overload of << operator
///@param s out stream
///@param d Distribution object to put in the out stream
std::ostream& operator<<(std::ostream &s, const Distribution& d)
{
	for(unsigned i = 0 ; i < d.size ; i++)
	{
		s << d.data[i] ;
		if(i < d.size -1)
			s << "\n";
	}
	return s;
}